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ABSTRACT 

We study the effects of an elastic crust on the oscillation spectrum of superfluid 
neutron stars. Within the two fluid formalism, we consider Newtonian stellar models 
that include the relevant constituents of a mature neutron stars. The core is formed 
by a mixture of superfluid neutrons and a conglomerate of charged particles, while the 
inner crust is described by a lattice of nuclei permeated by superfluid neutrons. We 
linearise the Poisson and the conservation equations of nonrotating superfluid stars 
and study the effects of elasticity, entrainment and composition stratification on the 
shear and acoustic modes. In both the core and the crust, the entrainment is derived 
from recent results for the nucleon effective mass. Solving the perturbation equations 
as an eigenvalue problem, we find that the presence of superfiuid neutrons in the crust 
and their large effective mass may have significant impact on the star's oscillation 
spectrum. 
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1 INTRODUCTION 

Neutron stars represent very complex systems, the understanding of which relies on uncertain theory combined with largely 
indirect observational evidence. The nature of the deep neutron star core, and the state of matter within, remains poorly 
constrained. However, it is clear that the quality of the gathered data, and our interpretation of it, is improving. As evidence 
of recent progress, with direct impact on the work discussed in this paper, we may consider the quasiperiodic oscill ations 
observed in the tails of giant magnetar flares ( Israel et al. 2005 : Strohmaver fc WattsI 20061 : Watts fc Strohmavei 20061 ) . The 
interpretation of the data in terms of torsional oscillations of the neutron star's crust provides us (at least in principle) 
with a handle on the associated physics. These observations have le d to a renewed interes t in neutron star seismology, 
and work addressing issues ranging from t he crust equation of state Jsteiner fc Watts | 200E ). to the role of the magnetic 



field (and th e coupling to the fluid core) ( Glampedakis et al. 20061 : van Hoven fc LevirJ I2OI1I : IColaiuda fc Kokkotasll201ll 



Gabler et al.ll2011) and even the relevance of the superfluid neutron s that permeate the nuclear lattice at densities beyond 
neutron drip (jSartiuelsson fc Andersson 20091 : Andersson et al. 20091 ). A similar breakthrough is hinted at by the observed 



cooling rate of the compact remnant in Cassiopeia A. The data provides clear evidence for the presence of superfluid neutrons 
in the star's core (associated with addi tional cooling du e to the Cooper pair b reaking mechanism), leading to a constraint on 



the superfluid transition temperature ( Page et al. 20 111 : Shternin et al. 201ll l. This constrains the theoretical models, which 



depend crucially on many-body interactions at extreme densities. Combined, the information gleaned about the crust region 
and the superfluid core provides us with a clear motivation to improve our seismology models accordingly. This is not a trivial 
task, but considerable progress has been made on the key aspects, realistic models would now seem to be within reach. 

The work reported in this paper takes a serious step to wards the modelling of r eal neutron star s eismo logy. Our effort 
should be considered in the context of previous work, such as McDermott et al. ( 19881 ) and Strohmavei (1991) who accounted 



* E-mail:andrea. passamonti@uni-tuebingen.de 



2 A. Passamonti & N. Andersson 



for the crust elasticity and finite temperature effect s, the modellins of superfluid dynamics by (for example) iLindblom fc Mendell 



1994 ) and most recently Passamonti et al. I 1I2OO9I) and Passamonti fc AnderssonI (|201ll ) in Newtonian gravity and 



Lin et al 



200s ) , Samuelsson fc AnderssonT i 200S ) and Kantor fc Gusakovl ( 2011 ) in the frame work of General Relativity, and finally the 



effort aimed at including the magnetic field (obviously a key aspect for magnetars!) ijColaiuda fc Kokkota I2OI1I : iGabler et al 



2OIII ). These studies provide key insights into the individual pieces of physics, but do not combine them, ft is obviously 



desirable to do so, and the present work takes important steps in this direction. We consider the oscillations of a neutron star 
model with an elastic crust permeated by superfluid neutrons and a fluid core with two dynamically distinct components. We 
develop, and implement, the relevant conditions at the interface between these two regions. We do not, however, account for 
the magnetic field. Neither do we (fully) account for thermal effects or consider the problem in General Relativity. These are 
obvious shortcomings of our final model, especially since we cannot meaningfully make use of a truly realistic equation of state 
unless we work in relativistic gravity. However, our model still represents the state-of-the-art for the coupled multi-fiuid-elastic 
aspects and the analysis provides valuable insights that will assist future developments in this problem area. 



2 THE TWO "FLUID" MODEL 

We want to build a neutron star model that accounts for the degrees of freedom associated with the crust elasticity and 
the superfiuid components, both in the crust region and in the core, ft is well-known that this involves modelling the core 
as a two-fiuid system, with the superfiuid neutron being treated as distinct from a charge-neutral conglomerate of protons 
and electrons. If the latter compon ents are strictly co-moving, and we ignore iss ues associated with superconductivity and 
the presence of magnetic fiuxtubes ( Glampedakis. Andersson &: SamuelssonI 201lh . then we may ignore magnetic field effects 
altogether. This is obviously not a true representation of a neutron star core, which should be magnetized, but the model is 
nevertheless useful as it allows us to explore the relevance of the additional degree of freedom implied by superfiuidity. 



2.1 The dynamics 

In the outer core of a neutron star, we expect superfiuid neutrons and superconducting protons to be present. Meanwhile, 
in the inner crust, at densities above the neutron drip threshold, a fraction of neutrons are superfiuid and coexist with the 
lattice of heavy nuclei. As we approach the neutron drip density, these "free" neutrons disappear and all the constituents of 
the crust are bound in nuclei. A two "fiuid" formalism is sufficient to describe these key aspect of neutron star dynamics, but 
we have to carefully consider the nature of each component in the crust and in the core. In the core it is natural to distinguish 
the superfiuid neutrons (denoted by n in the following) from a conglomerate of charged particles (protons/electrons, denoted 
by p). The distinction in the inner crust is not quite so obvious, see the discussion by , Andersson et a l. (2011). In this case, 
we consider the dynamics of baryons confined in the nuclei (denoted by c) permeated by a gas of free superfiuid neutrons 
(represented by f). 

With these definitions, the component fractions of the crust and the core are in general discontinuous at the crust/core 
interface, even if the total number of neutrons and protons is continuous. Therefore, the choice of junction conditions at the 
crust/core transition is very important for the oscillation dynamics of the star. 

The dynamics of a two-component star are governed by mass a nd moment um conservation equations for each constituent. 



^+V,(px^i) =0, (1) 



together with the Poisson equation for the gravitational potential (jPrixl 120041 ). The conservation equations take the form: 

dt 

{dt + viV,){v^ + e^w\^) + V.(m. + $) + e^wl^V.v^, = /f/px , (2) 

where and lOxy = vl^ — Vy are, respectively, the constituent and the relative velocities, while /ix = px/mx represents the 
chemical potential (we will assume that the neutron and proton masses are equal, i.e., we take m = nip = rrin). The mass 
densities are defined by px = rnn^, where nx is the constituent number density. The gravitational potential $ obeys the 
Poisson equation, 

V^'I' = 47rGp, (3) 

where p = px is the total mass density. In the interest of economy we will, whenever we write down an equation that 
describes both the dynamics in the core and the crust, imply that the two crust components have x either "c" or "f". In the 
core, these constituent indices should be replaced by "p" and "n" and the shear modulus should obviously be set to zero. 

To complete the dynamical equations, we need to provide the quantity Ex, that accounts for the entrainment, and the 
"external" force density /f that acts on each fluid component. As we will see later, the force can be used to account for the 
crust elasticity (and, in general, also the magnetic field). Finally, the system of equations is closed by an equation of state 
(EoS). It can be described by an energy functional that ensures Galilean invariance; 

£ ^£ {p^,'wly) . (4) 
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The chemical potential /ix and the entrainment parameter Ex are then determined from0 

(5) 



Px£x 



d£_ 



Py 

d£ 



dw'iy 



(6) 



2.2 Equilibrium model 

Let us first consider the equations that need to be solved in order to determine a non-rotating superfluid neutron star in 
dynamical and chemical equilibrium. We consider a static, nonrotating, background star in which matter is in /3-equilibrium, 
which mecins that jl = /ip — /^n- Moreover, we assume that the crust is unstrained, which means that the elasticity affects 
only the linear perturbations of the system. Under these assumptions, the star is described by the following equations: 

V/i = -V$, (7) 
= 47rGp. (8) 

Equation (O is obtained from the stationary Euler equations, the /3-equilibrium condition and the standard relation between 
the pressure and chemical potential VP = pV/i. The system of equations Q and ((SJ can be written as a single ordinary 
differential equation 

V^/i = -47rGp. (9) 

The background configuration is therefore completely determined once we choose an EoS that provides a second relation 
between the mass density and the chemical potential. More details on the EoSs used in this work are given below, in Sec. |3l 

2.3 Lagrangian perturbations 

In order to model the dynamics of the crust, where the elastic restoring force depends on the deviation from the unstrained 
configuration, it is natural to turn to Lagrangian perturbation theory. In the problem considered in this paper, the relevant 
perturbations are trivially related to the Eulerian variations of the various quantities, but in the general (rotating and possibly 
pre-strained) case th e Lagrangian descripti on becomes more involved. However, there has been recent progress in developing 



the required theory IjAndersson et al.ll201ll ). and we draw on these results in the following. In fact, our analysis provides the 



first application of the new formalism. 



In the two-fluid situation, it is natural to introduce a Lagrangian displacement vector for each fluid component (jAndersson. Comer fc Gr 



20041 ) . In order to describe the motion associated with these vectors, we define the Lagrangian perturbation of a general (scalar 



or vector) quantity Q along the vector field 

AxQ = 5Q + £5,Q, (10) 

where SQ is the Eulerian perturbation of Q and denotes the Lie derivative along . 

The dynamics of neutron star oscillations can be studied by linearising the system of equations HJ-®. The linearised 
mass conservation and Poisson equations are then given by; 

AxPx + PxV^Ci = , (11) 

V^^^ = 47rG ((5px + Spy) . (12) 

Meanwhile, for non-rotating stars, the linearised momentum conservation equations take the form ( Andersson et al. 20091 ) 

(1 - ^ + + ^^^^ + ^i^J^-^ - i^'^i) ^jf^- + V«Ac/ic = yV'a.j , (13) 

(1 - ^t) S + iS + + Vj V.cE. - (V,C0 V./if + V.Af/lf = , (14) 



where the elastic stress tensor is defined as 

2 



= A [y.^j + V,C- ) - -p (V^-^) 5., . (15) 
The shear modulus is denoted by /i, and must not be confused with the chemical potential ((5]). For nonrotating configurations 



^ In lPassamonti et al ] l|2009l ). the definition of the entrainment parameter ex given in equation (7) contains a typo. Here, we provide in 
equation the correct definition. 
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in /^-equilibrium, equations (|13|l and (|14|) assume the simpler form; 

As noted earlier, Equations (jlGp and p7|) also describe the dynamics of the neutron star core provided we replace, respectively, 
the indices c and f with p and n, and set jl to zero. 



2.4 Chemical gauge and crust/core junction conditions 



In the inner neutron star crust, the protons and a sizeable fraction of the neutrons are bound by the strong interaction. 
They form a lattice of nuclei that become heavier towards the crust/core transition, and which supports elastic stresses. The 
remaining fraction of (unbound) neutrons form a superfluid that may flow through the confined nucleons. This superfluid 
component first appears at the neutron drip density, pnd = 4.3 x 10^^ g cm""', which marks the transition between the inner 
and outer crust, where there are no free neutrons. 

In the inner crust there are, at least, two dynamical degrees of freedom; a confined component of protons and neutrons 
and the free superfluid neutrons. However, the distinction between superfluid and confined neutrons is not obvious. It depends 
on the local density and th e dynamical time-scale of the process under consideration. Adopting the operational definition due 
to Chamel fc Carter a neutron may be considered free when (on the time-scale considered) its energy is sufficient to 

overcome the potential barriers separating the nuclei, either classically or by quantum tunnelling. At the bottom of the crust, 
where the nuclei become denser and form exotic structure (the so-called pasta phase), the potential wells become closer and 
marginally bound states may exist. These states may penetra te the potential barriers on time-scales that are macroscopically 
long (but cosmologically short), see IChamel fc Cartel! |iooi) for discussion. 

Within the two-fluid forma l ism, this uncertainty may be expressed in terms of the chemical gauge parameter Uc introduced 
by jCarter fc Samuelssonlbood : IChamel fc Carterlbood : Undersson et al.lboill '): 



ni = Un + (1 - flc) rip, ric = acUp , (18) 

where nt and ric are, respectively, the number density of the free neutrons and the confined component, while rin and rip are 
the total number densities of neutrons and protons. On e can show that the neutron conjugate momentum and the chemical 
poten tial are independent of the chemical gauge choice ( Carter fc Samuelsson 20061 : Chamel fc Carter 20061 : Andersson et al 



2011f) as long as Oc is either held fixed or depends only on the nuclear charge number Z. This gauge independence is important 
for the derivation of the crust/core junction conditions. 

In order to study the oscillations of a given neutron star model, we must prescribe boundary conditions at the star's 
origin and surface. In a model with a crust, we must also specify conditions at the crust/core interfa ce and at the outer /inne r 
crust transition. In this work, we implement the static limit of the crust/core conditions derived bv I Andersson et al.l (|201ll ). 
From the mass conservation equations, we obtain the conditions for the radial component of the Lagrangian displacements at 
r = Rcc- 

= Cp, (19) 

cF = + (20) 

I — Xc 1 — Xc 

As in the single fiuid case the other junction conditions may be derived from the total momentum conservation equation, 



which leads to the following vertical and horizontal projections ({Andersson et al.ll2011 



(AeP+|/iV.Cc-2/i^)=0, 



N'N') (/iiV^-(V,C^+V,C,^)) =0. 



(21) 
(22) 



where Qij is the (fiat) metric. A*'' is a urut vector orthogonal to the crust/core interface and (...) denotes the change of 
a physical quantity across the interface. Equations (|21|l and (|22p reduce to the standard traction conditions in the single 
fiuid case. Finally, from the chemical gauge independence we can determine the condition for the chemical potential of the 
superfiuid neutron components; 

Ac/if = ApAn . (23) 
Prom the junction condition (|19p . it is clear that for nonrotating models in /3-equilibrium equation (|23p is equivalent to 

(5/Xf = Sfln- 

These conditions complete the description of the problem that needs to be solved in order to determine that star's 
oscillation spectrum (see Sec. 14.21 for the inner/outer crust boundary conditions). 
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3 MODEL EQUATION OF STATE 

In order to make maximal use of the theoretical framework, we ought to build our models using a realistic equation of state 
for supranuclear matter. We will, however, resort to using a simple model EoS. There are several reasons for this. Most 
importantly, we do not yet have a realistic EoS that provides all the different parameters required in our analysis. Most 
tabulated EoS simply provide the pressure vs density relation, whereas we need (at the very least) detailed information about 
the composition, the sup e rfluid pairing gaps and the entrainment between neutrons and protons. Progress in this direction is 
being made, c.f., Chamel ( 20081 ) . but the models are not yet complete. Still, the developments have reached the point where it 



would be meaningful to consider a fully relativistic analysis (as is needed to make the use of a realistic EoS meaningful). This 
is, however, beyond the scope of the present work. Our main interest here is to probe the key issues and consider the associated 
phenomenology. For this exercise it makes sense to focus on a simple analytic model with freely adjustable parameters. 

We model the neutron star matter in terms of two polytropic EoSs. Despite the simplicity of the construction, combinations 
of polytropes enable us to study the effects of composition stratification and symmetry energy on the oscilla tion spectrum. 
Moreover, we have already used these EoSs in previous work ( Passamonti &: Anderssonlboil : Siderv et al.ll201ol ). which means 



that we can build on, and compare to, the previous results. As a new application of the polytropic EoS, we construct stratified 
two-fluid models with different parameters in the crust and the core in order to approximate the qualitatively different 
compositions in these two regions. 

In most astrophysical systems the relative velocity between the two fluids is small. Therefore, equation (j4| can be expanded 
in a series: 

£ = £o {pi, Pc) + ao {pi, pc) lUfc + O {wt) , (24) 

This approximation has the advantage that the bulk EoS, S^, and the entrainment parameter, qq, can be independently 
specified at w}^ — 0. From equation ((SJ it follows that the entrainment parameter £x is related to the function qq by 

PxEx = 2ao . (25) 

In a co-moving background, the density perturbations can be expressed in terms of the chemical potential perturbations as: 

Ape = SaApi + SicApc , (26) 
Apc = ScfApt + 5ccA/ic , (27) 

where 



3.1 The bulk equation of state 



Gi ven the above re l ations, we need to f i rst of all provide the energy S p. The first EoS we consider has also been used 
bv lPrix et J] l|2002h : IVoshida fc Eriguchil (|2004h : IPassamonti et al.l !\2QQ^ ): 



So 



K 



2Ka, 



-L (1 ~1~ ^ syin) Xc 1 (1 O sytn) Xc Xc [1 (1 O syva) Xc\ 

where K is a. pol ytropic constant, Xc is the confined nucleon fraction and (Jsym is a param eter that can be related to the 
symmetry energy ( Prix et al. 2002). In this EoS, both Xc and asym are taken to be constant ( Passamonti et al. 20091 ) . Using 
equation (|28|) and (|29p . we readily obtain 



-Pipe + 



(1 + 2(7, 



J ^2 
-pc , 



(29) 



5ff 
<5cf 



2K{1 + asym) 

1 -|- CT sym " (1 H" 2(Jsym) 



(30) 



(31) 
(32) 



2K (1 -I- asym) 

Sfc — asym *5cc . 

As discussed by Passamonti et al. 1 20091 ). this EoS leads to non-stratified stellar models. The global properties of these models 
are similar to those of the A'^ = 1 polytrope in the single fluid case. For instance, the pressure is related to the total mass 
density by the standard equation P — KfP" . The set of models constructed from this EoS will be referred to as models A, c.f. 
Passamonti et al. ( 20091 ) . 

At this point, it is useful to pause and make some comments on the numerical implementation and the resultant stellar 
models. In our numerical code, we use dimensionless units based on the gravitational constant G, the central mass density po 
and the stellar radius R. Unless directly specified, all models A considered in this paper are such that the crust /core transition 



is at _Rec = 0.9ii. Therefore, from the expected crust/core transition density pcc = 1.2845 x 10 g cm (jDouchin fc Haense! 



200 ll ). we can determine the central mass density of the star. At r 
g cm 



Kcc the dimensionless mass density is pcc = 0.109 po, 
hence po — 1.175 x 10^^ g cm~"^. Meanwhile, the dimensionless mass of models A is Mj {poR'^) — 1.273. Therefore, the physical 
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Figure 1. This figure shows some properties of a two fiuid polytropic model D with Xc = 0.3 at r = Rcc- The left panel displays the 
radial profile of the mass density and the chemical potential. The right panel shows (on a semi-logarithmic scale) the fraction of the 
crust /core constituents (vertical axis) and the total mass density (horizontal axis). In the right panel the vertical dotted lines denotes 
the density of the crust/core transition pcc and the neutron drip density pnd. Note that the total number of neutrons and protons is 
continuous at r = Rcc, as in the crust a fraction of neutrons is confined in the nuclei of the lattice. 



Table 1. This table provides the main parameters for the stellar model sequence D. The core is described by the EoS 1 133 II with A'n = 1, 
A^p = 1.4 and fcn = 0.766, fcp = 2.256. The crust/core transition has been set to Rcc = 0.88iJ. In the crust, we choose the same polytropic 
indices for the free neutrons and confined nucleons, respectively, Nf = 2.9 and Nc = 1.2, but vary at r = Rcc- For different values 
of Xc at the crust /core interface (first column), we show in the second and third column the dimensionless EoS parameters kc and fcf, 
respectively. The dimensionless mass is given in the fourth column. Note that the coefficients fcx are given in units of GR^p^ 



Xc 


kc 


k( 


M/ (poR3) 


0.1 


3.574 


0.258 


1.149 


0.3 


1.431 


0.281 


1.155 


0.5 


0.935 


0.316 


1.161 


0.7 


0.706 


0.377 


1.167 


1.0 


0.525 


13.172 


1.176 



mass and the radius can also be determined. For instance, in the case of the canonical neutron star mass, M = 1.4 Mq, we 
obtain the radius R = 12.29 km. 

It is obviously not the case that the component fractions are constant in a real neutron star. We can account for 
composition gradients by considering another simple combination of polytropes; 

So = fcfp^f + fccp? , (33) 

where fcx and 7x are constants. This EoS does not have the symmetry energy term, but it can leads to models with varying 
composition when 7f 7^ 7c. From equations ((5]) and H33|) it follows that the chemical potential and the corresponding mass 
density are related by 

Px-(^) , (34) 



. fcx7x , 

where the polytropic index is given by A^x = (7x — From this result we can determine the proton fraction for a given 

stellar model by imposing /3-equilibrium. After some calculations, we obtain: 



{jckc 



-JVf 



(35) 



For /i — >■ 0, it is clear from equation (|35|l that Xc vanishes for A^f > A^c and tends to unity when A'f < Nc. 

In a "realistic" neutron star, the proton fraction decreases from the centre to the crust/core interface, while the confined 
component increases from the bottom of the crust and reaches unity at the neutron drip. This behaviour can be approximated 
using the EoS (|33[) provided we choose different values for the polytropic indices in the core and the crust, i.e. let {Nn, Np) 7^ 
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Figure 2. This figure shows the dependence of the effective masses on the baryon number density. In the left panel, we show the neutron 
effective mass for both the co re (neutrons) and the crust (free neutro ns). In the crust, the filled circles denote the data set determined 
by Chamel and collaborators llCarter et al.ll2005l : IChamell 1200^. 12006*). These data have been fitted with the curve (dashed-line) shown 
in the left panel. For model D with Xc{Rcc) = 0.3, we display in the right panel the effective mass for the charged component (core) and 
the confined particles (crust). 



{N{,Nc). For simplicity, we also assume that the equilibrium star is such that the total mass density, the chemical potential 
and its first derivative are continuous at the crust/core interface. We have found a good compromise that satisfies all these 
conditions by choosing A'n = 1, A^p ~ 1.4 and fcn ~ 0.766, kp — 2.256, in the core and A^f = 2.9, A^c = 1.2 and kt — 0.281, 
fee = 1.431 in the crust. The crust/core transition has been fixed at R^c ~ 0.88-R. Note that the coefficients fcx are given in 
units of GR^Pq ''^ . With these parameters, the proton fraction in the core ranges from Xp = 0.2 at r = to a;p = 0.1 at 
the crust/core boundary. These values are on the high side, but this is not too concerning here. In the crust, the charged 
component varies from Xc = 0.3 at r = Rcc to a;c = 1 at the inner/outer crust interface, r = i?ND- In this model, the total 
mass density at the interface is pcc = 0.115 po in dimensionless units. As in the case of model A, we can therefore determine 
the central mass density, po = 1.117 x 10^^ g cm~^. For a star with mass M = 1.4 Mq, we can determine the stellar radius 



from the dimensionless mass of this model; from M/{poR^) = 1.155 we find that the radius is J? = 12.92 km. 

Fig. H] illustrates the properties of this crust-core model. The mass density and the chemical potential are shown in the 
left panel and the constituent fractions are provided in the right panel. The crust /core and neutron drip transitions are also 
indic ated. We have assumed a neutron drip density pnd = 4.3 x 10^^ g cm~^ I Shapiro fc Teukolskv 19831 : Douchin fc Haense] 
2001f ) which for this stellar model corresponds to TiND = 0.99887?. For the same core model, we have also studied solutions 
with different values of Xc at the crust/core interface. From equation (|35|l and the continuity conditions for the total mass 
density and chemical potential, it is evident that the coefficients fcc and kf depend on the parameter a;c at r = Rcc- Their 
values and the total mass of the stellar model are listed in Table [T] We refer to t his sequence of stratified stars as models D 
in order to distinguish them fr om models B and C considered in previous work ( Passamonti et al. 20091 : Siderv et al. 20ld : 
Passamonti &: Andersson 2011). 



3.2 Entrainment 

In a real neutron star, superfluid neutrons may flow through the lattice of nuclei in the crust and also coexist with supercon- 
ducting protons in the core. The different constituents are coupled through entrainment, mutual friction, crust/core boundary 
conditions, EoS, etc. The entrainment is a non-dissipative process that arises from strong interaction. Owing to this effect 
each superfluid constituent can carry along a fraction of the other component. From a dynamical point of view, it is natural 
to account for the entrainment by introducing an effective mass associated with each component of the system and which is 
in general different from the bare particle mass. Meanwhile, in the two-fluid formalism, the entraiment is usually described 
by the parameter Ex defined by equation ([6]). It is related to the effective mass mj as 

Ex = 1 - ^ . (36) 
m 

Entrainment affects the constituent conjugate momentum in such a way that it is not longer aligned with the transport 
velocity of the fluid constituent. 
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In recent work, IChamel l|2008l ) determined the nucleon effective masses in a neutron star core. The results are based 
on the two-fluid formalism and describe the strong interaction between nucleons using a two-body force of the Skyrme type 
with diff erent paramet risations. In a non-relativistic two-fluid model the effective masses can be described by the following 
relation (jchamellbooj ): 



1 



" 1 + P-iph 



where ph is the baryon mass density, and /Sa is defined by 



3 k 2B3 
m 



(37) 



(38) 



where is the reduc ed Planck's constant. The quantity Bs is a coefficient that depends on the parametrisation of the Skyrme 
force (jChamellbooi l. 

In this work, the superfluid neutron star core is modeled by two fluid polytropes and the effective neutron mass is 
determined by equation (I37|l . For the proton effective mass, we use the EoS constraint (|25|l . i.e. 

Pn£n — ppEp , (39) 

together with the definition psp . This model is obviously not consistent, as the effective mass relation (|37|l has been determined 
for "realistic" EoS and not for polytropes. However, this is the first time that the complete oscillation spectrum of superfluid 
neutron star with an elastic crust is considered. Hence, it is natural to focus on the phenomenology sacriflcing (to some extent) 
realism. Our stellar model has the key features expected of a real neutron star, but the "proportions" may not be quite realistic. 
This should, however, be easy to fix once the problem is considered in General Relativity. The relevant developments are, in 
fact, already under way. 

The effective mass of the free superfluid neutrons in the inner crust has been determined with in the band theory approx- 
imation, but only for a few speciflc densities ( Carter. Chamel fc Haensell 2005 : ChameJ 20051 . 2006). The results indicate that, 
in the deep crust, the effective mass can be as large as rrif = 15 Arm. As already discussed bv I Andersson et al. (200^), who 
studied the problem within a local (plane-wave) approximation, these large entrainment values may have considerable effect 
on the properties of the oscillation spectrum. The currently available entrainment data is shown in the left panel of Fig. (2] 
We interpolate these data points with a curve, and determine the proton effective mass with equation (|39p by applying the 
method that was used for the core. The results are shown in Fig. [21 where the right panel displays the proton effective mass 
for a model D with Xc = 0.3 at the crust/core interface. 



4 THE EIGENVALUE PROBLEM 

For non-rotating stellar models, the study of the linearised equations (|f f p - (|I2p and (|f 6p - (|I7p is considerably simplified if we 
expand the perturbation variables in vector harmonics. For each set of harmonic indices {I, m) (associated with the familiar 
spherical harmonics V,™), the relevant equations then depend only on the radial coordinate r and time t. In order to set up an 
eigenvalue problem, we assume that the perturbation variables have a harmonic time dependence e"^' (this essentially amounts 
to taking a Fourier transform of the equations). The mode eigenfrequency a is in general a complex number, where the real 
part describes the oscillation frequency and the imaginary part the growth/damping time of the mode under consideration. 
However, in this work the eigenfrequency a is real since we are not accounting for dissipative processes or mechanisms that 
may trigger instabilities. 

The perturbations can be decomposed into two classes of different parity, namely the polar (spheroidal) and the axial 
(torsional) perturbations. The spectrum of non-rotating stars is also degenerate with respect to the azimuthal harmonic 
m. This degeneracy is removed when we account for additional physics like rotation or the star's magnetic field. Oscillation 
modes of a non-rotating star are commonly classi fied by the dom inant restoring force that acts on the displaced fluid elements. 



the harmonic index I and the radial number n l|Cowling||l94lh that labels the "overtones" of a given set of oscillations. In 



simple models the latter parameter can be associated with the number of radial nodes in the relevant eigenfunctions, but this 
identiflcation no longer holds for more complex systems. 

The axial part of the oscillation spectrum is comparatively simple. For non-rotating stars without magnetic fields and 
crust, the axial perturbations form a degenerate set of zero frequency modes. The crust elasticity breaks this degeneracy 
as it sustains a set of torsional shear modes. We denote the torsional modes (t-modes) by the symbol 't„, where I is the 
harmonic index and n labels the overtones. For any I, the frequency of the torsional modes assumes its lowest value for the 
fun damental torsional mod e and increases with n. In compact stars, the typical frequency of the ^to-mode is a few tens of 



Hz (McDermott et al. 198: 



The polar oscillation spectrum is more complex. In addition to the shear waves, the presence of pressure and com- 
position/thermal gradients may generate acoustic and gravity waves. The acoustic modes are mainly restored by pressure 
variations (p-modes) and cover the high frequency band of the spectrum, above I kHz. At lower frequencies, typically below 
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100 Hz, composition and thermal gradients can sustain the class of gravity modes (g- modes), which are restored by buoyancy. 
For a given I, there is an infinite series of pressure (gravity) modes with n ^ 1, whose frequency is increasing (decreasing) with 
n. The i,^ eigenfunction of the pressure, 'p„, and gravity modes, 'g„, have n radial nodes. The fundamental mode (f-mode), 
'f, separates the class of p- and g-modes. Its frequency scales with the averaged stellar density and no radial nodes are present 
in its eigenfunction. 

The crust elasticity affects the polar oscillation spectrum in two ways. First ly, it generates a class o f spheroidal shear 
modes (s-modes) and secondly it slightly modifies the acoustic mode frequencies ( McDermott et al. 19881 '). Spheroidal shear 
modes are represented by 's„, with n ^ 1. Their frequencies increase with n, and for small n they lie below the f-mode 
frequency. In the very low frequency band of the spectrum, there are also interface modes,which are denoted by 'i. These 
modes arise from the presence of internal interfaces and they strongly depend on the physical conditions near these interfaces. 
For any I, there is an 'i mode associated with each internal interface. 

The mode classification presented above allows us to describe the complete oscillation spectrum of non-rotating, non- 
magnetised single-fluid stars. In superfluid neutron stars, the multi-component dynamics of the system lead to additional 
degrees of freedom and the oscillation spectrum is therefore richer. In fact, in a two-fluid model we can (usually) decompose 
the dynamics into a co-moving (counter-moving) motion, where the fluid elements oscillate in (and out of) phase. As a result, 
we can distinguish two classes of acoustic modes that are commonly called "ordinary" and "superfluid" modes. The ordinary 
modes have spectral properties similar to the acoustic modes of single fluid stars, as their dynamics is mainly dominated by 
the CO- moving degrees of freedom (and they are therefore restored by the total pressure gradients etcetera). Meanwhile, the 
superfluid modes exist only in multi-component stars, and the counter-phase motion of the perturbed fluid elements dominates 
their dynamics. Apart from in particular cases, a neat decomposition of co- and counter-moving motion into "ordinary" and 
"superfluid" modes is not possible. A typical acoustic mode tends to exhibit both degrees of freedom. Nevertheless, in order 
to distinguish between these two families of acoustic modes, we denote the ordinary (superfluid) f- and p-modes with the 
symbols 'f° and 'p° {'■P and 'p^), respectively. 

The spectrum of superfluid neutron stars exhibits another important difference with respect to the standard single 
fluid mode ls. Superfluid neutron stars do not support an independent class of g-modes, even if the model has composition 
gradients 1 19951 : lAndersson fc Comeil I2OO1I : Iprix fc RieutordI I2OO2I ) . In a perturbed two-fluid system, deviations from 
chemical equilibrium tend to excite the counter -moving degree of freedom . The buoyancy only introduces small corrections 
to the superfluid acoustic mode frequencies (see lAndersson fc Comeill200ll . for more details). 



4.1 The polar perturbation problem 

The expansion of the polar perturbations in vector harmonics is given by the following expressions: 

Q^W''ir)Yim, & = V''{r)deYi^, = ' 5$ ^ 5^r)Yi^ . (40) 

All scalar perturbations are expanded in terms of the spherical harmonics Yim, c.f., the perturbed gravitational potential. The 
problem becomes more tractable if we work with dimensionless quantities. In the core, we choose to work with the following 
set of perturbation variables: 

Sflp W 5 fin 5$ 5$' 

yi = y2 = y-i = y4 = ys = — ye = — (41) 

r rg r fg 9 

where g is the gravitational acceleration. For spherical stars, this quantity can be related to the chemical potential by the 
equilibrium equation g = ^' — —jl'. 

The core is described by a system of six perturbation equations that we can write (formally) as 

^ = n[y.,a^6(r)] i,fc = 1...6, (42) 

where the variables yi are defined in pif) . and the function b{r) denotes a generic background field. The complete system of 
linearised equations is given in Appendix lAll 

In the crust, the perturbation equations are evidently more complex due to the presence of the elastic strain tensor. 
As in single fiuid models, it is important to choose a set of perturbation variables that makes the solution of the problem, 
including the implementation of the crust/core junction conditions, easier. In the inner crust, where a mixture of superfluid 
free neutrons coexist with a lattice of heavy nuclei, we deflne the following set of eight dimensionless perturbation variables: 

Zl = Z2 = PcA^c + o A*Vi4c - ^fJ-W^ zz = — 2:4 = — K 1 

r Qc \ i ) r gc \ r r 

_ W' _ Aflf _ 5$ _ 5$' 

Z5 = Zfi = zr = — zs = (43) 

r rg rg g 

where gc has the same units as the shear modulus, and has been introduced to make Z2 and 24 dimensionless. In this work. 
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we choose qc = Po, where po is the pressure at the origin r = 0. In the function Z2, the divergence can be expressed in terms 
of the harmonic expansion variables since 

V.Cc = W^ + -W,~ K . (44) 

r r 

The variable Z2 and Z4 are directly related to the radial and angular traction components (see Sec 12.41 ). This simplifies the 
implementation of the junction conditions. Another difference with respect to the functions used in the core {yt) is that Z2 and 
26 depend on the Lagrangian perturbation of the chemical potential. This choice provides a simpler boundary- value problem 
and less complex equa tions in the core. In the single fluid limit, i.e. when Xc = 1, the perturbation variables (|43|) reduce to 
the quantities used by McDermott et al. ( 19881 ) provided we set = pH Naturally, the variables 25 and za are obsolete in 
this limit. 

The linearised equations of the crust form a system of eight ordinary differential equations that we can write (again, 
formally) as: 

^ = Zk[z„a\fi,b{r)] i,fe = 1...8. (45) 
ar 



We provide the complete form of these equations in Appendix IA2I 

In the outer crust, all the neutrons are confined to nuclei, i.e. we have X{ = 0. Therefore, the perturbation problem reduces 
to the single-fluid case. The oscillation dynamics is then described by a system of six ordinary differential equations, where 
the variables 25 and Z(, are no longer used. 



4.2 Boundary Conditions 

To solve equations (|42p and (|45p as an eigenvalue problem, we must impose regularity conditions at the centre, boundary 
conditions at the surface of the star, and junction conditions at the two internal interfaces. The first internal interface marks 
the transition between the outer and the inner crust, and the second represents the point where the inner crust gives way to 
the fluid core. 

At the origin (r — 0), the regularity conditions lead to the following relations: 

I J/2 = <^ yi-\ — cr ys - 1/6 , (46) 

go go 

hi = — o- yi H o- ys-ya, (47) 

go go 

lyl = ye, (48) 

where go and i/° are, respectively, the gravitational acceleration and the perturbation variables at the centre. 

In Sec. 12.41 we determined the junction conditions at the crust/core surface r = Rcc- In terms of the dimensionless 
perturbation variables these conditions take the form 

^1 = yi, (49) 

22 = —rg [x^yi - Xpy2 + {xp - Xc) yi] , (50) 
qc 

24 = 0, (51) 

Xp Xc 1 ~ Xp /ro\ 

Z5 = -f yi + -. -ys. , (52 

I - Xc 1 - Xc 

z(i = yi - Z5, (53) 

Z7 = ys, (54) 

zs = ye ■ (55) 

The boundary conditions that must be satisfied at the inner/outer crust interface are the continuity of the radial dis- 
placement 2i, the traction components 22 and 24, the gravitational potential perturbation 27 and its radial derivative zs- 
The component of free superfluid neutrons vanishes on this surface. Given this, we impose a zero boundary condition for the 
Lagrangian perturbation of the free neutron chemical potential (A/if), i.e. take 25 — 0. 

Beyond the transition to the outer crust, the problem is equivalent to the single fluid approximation. Therefore, at the 
stellar surface we must impose the continuity of the traction perturbation and a co ndition for th e gravitati onal perturbation . 
The r — >■ -R limit of the linearised Poisson equation leads to the following expression (jUnno et al ] [l989, : ,Prix fc Rieutordll2002l ') : 



+ ^ri^* = -r lim p{r) , (56) 

R 1 — >R- 



In lMcDermott et al.l l|l98^ , the definition of the variable 22 contains a typo. The total mass density p must be replaced by the pressure 
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For a model such that the density vanishes at the surface, this means that we should have 

5$' + ^-^5$ = . (57) 
H 

At the star's surface, the perturbation variables must then satisfy the following relations: 

zi = 1 (58) 

Z2 = (59) 

Z4 = (60) 

Z8 = -(1 + 1)Z7, (61) 

where we have imposed a normalization condition for the confined nucleon Lagrangian displacement, zi = 1. 

The outer crust of the stellar models studied in this work is very thin. Neutrons start to drip out from the nuclei at 
r — Rnd — 0.9988-R. The mode frequency calculation is not really affected if we impose the boundary condition A/if — 
(zfj = 0) relevant to the inner/outer crust interface at the surface of the star. This is not surprising. The fact that the low 
density region has little effect on the global oscillations of the star also means that we can ignore the thin ocean that should 
exist below the density at which the crust is expected to melt (at finite temperatures). 



4.3 The axial perturbation problem 

Axial perturbations of spherical stars describe torsional oscillations which are sustained by the crust's shear. The corresponding 
vector harmonic expansion of the Lagrangian displacement is given by 



sin 6 



^^ = -U''{r)deYi„ 



(62) 



while all Eulerian perturbations of the scalar variables are zero. From equation (|62p . it is clear that the vector has zero 
divergence, Vi^x = 0. Therefore, we have Apx ~ A/ix = 0. 

Substitution of the expansion (|62|1 into the linearised equations (|16[1 -(|17 [) leads to a system of two equations: 



2 - 



Pc 



U'' 



r 



Ui = 

1 - Ef 

where we have defined 



Uc 



1 - En 



(63) 
(64) 

(65) 



As before, we have assumed a harmonic time dependence with frequency a, as before. The system of equations (I63|l - (|64p 
describes the linear dynamics of torsional oscillations in the inner crust. In the outer crust, there are no free superfluid 
neutrons, which means that we have Xf — 0. Therefore, equation (|64p becomes obsolete, while equation (|63|l reduces to the 
single fluid approximation provided Xc — 1 and the entrainment is set to zero, i.e. £x = 0. 

It is worth noting that, in both the outer and inner crust, equation (|63|l depends only on the Uc variable, which therefore 
completely deter mines the torsional oscillation m odes. That this is the case also in General Relativity has already been 
demonstrated by Samuelsson fc Andersson 1 2009f ). Once the eigenvalue problem is solved for the confined nucleons, we can 
determine the eigenfunctions for the free neutrons in the inner crust from equation (|64p . 

As in the polar case, we introduce a new system of dimensionless variables in order to simplify the boundary value 
problem. These variables are defined as: 



Sl 



S2 



<3c V r 



(66) 



where S2 is proportional to the transverse component of the traction. Equation (|63p now becomes a system of two first order 
ODE: 



dsi 
dr 

dS2 

dr 



= — S2 , 
P 



pc 2 2 

- — r a 

£* 



£1 

qc 



3S2 , 



(67) 
(68) 



The relevant boundary conditions are; the continuity of the traction at the crust/core interface, inner/outer crust interface 
and at the stellar surface. We therefore impose S2 = both at the crust/core boundary and at the star's surface. In addition, 
we normalize the radial displacement by setting si = 1 at the stellar surface. 
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Figure 3. For non-stratified models A with zero entrainment, this figure displays the 1 = 2 interface mode and the first five I = 2 
spheroidal shear modes. In both panels, the numerical results are shown as lines, while the values estimated from the plane-wave 
equation I I70I I are represented by symbols. In the left panel, we show the dependence of the interface and shear mode frequencies (vertical 
axis) on the elastic modulus parameter /i« (horizontal axis) for a model A with constant Xp = 0.1 and Xc = 0.3. In the right panel, 
we study the dependence of the spectrum on the confined component fraction Xc (horizontal axis). In this case, the models A have 
fit = 10~^. In both cases, equation l|70| l describes very well the spectrum of the spheroidal shear modes. 



5 RESULTS 

Let us now consider the detailed oscillations of the multi-component neutron star model developed in the previous sections. 
We focus on non-rotating, polytropic two-fluid models, and study the effects of varying component fractions, crust elasticity 
and entrainment on both the shear modes and acoustic waves. Building on the discussion in Sec.|4l we decompose the problem 
into polar and axial perturbations. The spectrum of each family of perturbations can be studied as a boundary value problem 
by solving the system of equations (|42[l and (|45p for the polar modes and (|67p - (|68p for the axial modes. In both cases we use 
a shooting method to a fixed point. For the polar sector, we set the fixed point at the crust/core interface, where the core 
and crust solutions are matched using the junctions conditions discussed in Sec. 12.41 and l4. 21 



5.1 Shear modes 



In simple plane-wave models, the frequency of the crust shear waves is proportional to the shear velocity; 



(69) 



However, in the case of a two- fluid model the shea r velocity depends on t he mass density of the confined nucleons and not on 
the total mass density, as in the single fluid case ( Andersson et al. 20091 ). Still, considering typical parameters, it is easy to 
see that the shear waves have low frequencies. One would also expect the shear model to be predominantly transverse. They 
should be present in both the polar and the axial sector, and we will refer to the two classes as spheroidal and torsional shear 
modes, respectively. The crust elasticity will also affect the acoustic mode frequencies by introducing very small corrections 
(roughly of the order of the ratio between the shear wave speed and the speed of sound) to the acoustic modes. Furthermore, 
the free superfluid neutrons in the inner crust may interact with the lattice of nuclei through entrainment and modify the 
spectrum. The main impact of the entrainment is expected in the deep crust, where the effective mass of the free neutrons 
may be very large (jChamcl 2005), c.f.. Figure [2] 

From a plane-wave analysis of the two-ffuid proble m, we can derive the f requency dependence of the shear modes on the 
mass density, the entrainment and the shear modulus ( Andersson et al. 20091 ): 



(70) 



where k is the magnitude of the wave-vector. We are interested in testing the accuracy of this rough approximation against 

full numerical results for stratified and non-stratified superflu id models. 

The shear modulus of the crust can be determined from i Strohmaver et al]|l99ll ): 

0.1194 m {Zef 



i-f 0.595 {To/ry 



(71) 
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Figure 4. For a model A with constant Xp = 0.1, Xc = 0.3, and /i* = 10~*, we sliow in this figure the frequency dependence of the 
I = 2 interface and spheroidal shear modes on the entrainment parameter e* (horizontal axis). The numerical and estimated frequencies 
are represented using the same notation as in Fig. \3\ The results show that the shear mode cigenfrequencies are well approximated by 
equation jTOj. 



This result was inferred from Monte Carlo simulations of the Coulomb interactions in the neutron star crust. The quantity Ui 
is the ion density, Z the atomic number of the nuclei and a — {3/4nni)^^''^ is the average inter-ion spacing. The ratio between 
the Coul omb and thermal energy is pa rameterised by F = (Ze)'^ /aksT, where T is the temperature and fcs the Boltzmann 



constant. 



Farouki fc Hamaguchil ( 1993[ ) have discussed plausible values for the parameter Fq. 



To be truly consistent with (|7ip . a stellar model should obviously be obtained from a realistic EoS. As we have already 
discussed, we adopt a different strategy in this work. We focus on simple polytropic two-f luid models and determine the shear 
modulus from a phenomenological relation. Detailed calculations for realistic EoSs (see iDouchin fc Haensell (|200in ) suggest 
that the specific shear modulus is almost constant across the crust, i.e. 



P 



10 cm s 



(72) 



Hence, it is reasonable to (as a first approximation) consider simple models such that 

/i = M*P, (73) 

where ^i, is a constant parameter. We consider only this simple model, and in the numerical code, we use its dimensionless 
form jit, which is defined as follows: 

Ai, = A* GpoR^ . (74) 
The constant parameter (|74|) may be estimated at the crust /core interface by using equation (|72p . the transition mass density 



Pec 



1.2845 X 10^* g cm at r = iZcc, and assuming typical neutron star values for the central mass density po 



10^^ gcm"^ 



and radius R — IQ km. This way we obtain /i* ~ 10 &t r — Rcc- 



5.1.1 Spheroidal shear and interface modes 

We consider first the spheroidal shear modes of the set of A models; stellar configurations with constant component fractions. 
In order to explore the dependence on the shear modulus, we focus on a model with Xp = 0.1, Xc = 0.3 and zero entrainment. 
The resulting frequencies of the I — 2 shear modes and the interface mode are shown in the left panel of Fig. [3] for different 
values of the shear modulus (/t,). We compare our numerical results to the plane- wave relation (|70p . To do this we first 
"infer" the effective value of the wave vector k in equation 1)70^ for a single numerical data point {a, jl) and then use this 
value in equation (|70|) to extend the relation to other parameter values. From the results in Fig. [3l we see that the numerical 
results scale according to the phenomenological relation. The agreement is, in fact, very good for the shear modes. Meanwhile, 
the interface mode frequencies deviate from the plane-wave relation. This is expected since these modes have a different 
origin; they depend sensitively on the physical conditions near the crust/core interface. Interface modes typically have small 
frequencies and amplitudes concentrated around the internal interfaces. A characteristic property of these modes is a cusp 
in the radial displacement eigenfunction at the internal boundary layer. In this work, we do not consider neutron stars with 
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Figure 5. This figure displays the 1 = 2 spheroidal shear mode eigenfrequencies for model A (left panel) and D (right panel). In both 
cases, the shear modulus parameter is fit = 10~*. On the horizontal axis, we show the radial number n of the shear mode ^Sn, while 
the vertical axis displays the dimensionless mode frequencies. The frequencies for the zero-entrainment model are shown as filled circles, 
those corresponding to the Chamel entrainment as filled diamonds, while the open circles are the values estimated from equation l|76|l . 
We show also the shear mode frequencies (denoted by stars) for models without free superfluid neutrons in the crust {xc = 1). 



an ocean, which gene rates another interface mode associated with the ocean/crust transition (see McDermott et al. 198 



Piro Sz Bildstenll2005l . for more details). 



Let us now apply the same strategy to the dependence on the fraction of confined nucleons Xc- In this case, we consider 
a model A with Xp = 0.1, zero entrainment and /i» = 10~^. The results shown in the right panel of Fig. |3]confirm the good 
agreement between the numerical shear-mode frequencies and the analytical relation (|70p also in this case. 

We can also study the dependence on the entrainment parameter £*. To do this we construct Model A to be a star with 
Xc = 0.3, and shear modulus /i, = 10~*. We consider two distinct entrainment configurations. In the 

from the relation £x = e/xy . We have already considered this 



first case, we set the parameter e constant and determine 

configuration in a different context ( Passamonti et allboogl : Passamonti fc Andersson 



the entrainment model described in Sec. 13.21 which is based on the work of Ichame 



2011 



i; In the secon d case, we investigate 



I2OO5I . bood . I2OO8I ). When e is taken 



to 

be constant, equation H70p describes the frequencies of the shear modes very well, see figure U However, the plane-wave 
approximation is not accurate for the model with a realistic entrainment profile. This result is (obviously) expected, as 
the plane- wave approximation follows from a local analysis (based on constant parameters). Nevertheless, we find that an 
"average" of equation (|70p over the crust reproduces the shear frequencies of this model rather well. This means that we can 
use the following shear mode relation: 

2 _ {£*)o 



(pc ) av 



(75) 



where ( )av denotes an average over the crust. In the left panel of Fig. [SI we compare the shear mode frequencies of model A 
without entrainment, i.e. with e = 0, to a model based on "realistic" entrainment. The results show that the relation between 
these two sets of data is accurately described by equation (|75p . i.e., 

o- (e*) = V (e*)ai.cro , (76) 

where ao is the frequency for the zero entrainment model. 

Next we consider the effect of composition gradients on the shear modes, making use of the D models introduced in 
Sec. [3] With two distinct values for the polytropic indices of the core and the crust, these models mimic the composition 
of realistic EoS, i.e., the proton fraction decreases from a central value, while the fraction of confined nucleons in the crust 
increases away from the interface. We study a sequence of models with the same core configuration, such that Xp = 0.2 at 
r = 0, and different constituent fractions in the crust. More details on these models D are given in Sec. [3] Before quantifying 
the effect that a composition gradient has on the shear modes, we consider the realistic entrainment profile for these stratified 
models. To this end, we select a model D with Xc = 0.3 at the crust/core interface and let /i* = 10~*. The obtained shear 
mode eigenfrequencies are shown in the right panel of Fig.O Approximate frequencies determined from the formula (|76l) . are 
(again) in good agreement with the numerical values. 

In order to emphasize the importance of entrainment on the spectrum, we also show in Fig. [5] the shear mode frequencies 
for stellar models without free superfluid neutrons in the crust (xc = 1). With a large effective mass, the relative motion 
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Figure 6. This figure shows the efli'ects of composition gradients on the spheroidal shear modes. The background stars are described by 
models D with different Xc at the crust/core interface, r = Rcc, but with the same shear modulus, fit = 10~*, and zero entrainment. In 
the left panel, we show the dimensionless eigenfrequencies of the first six I = 2 spheroidal shear modes (lines) and their values estimated 
with equation l|78|l (symbols). The accuracy of this equation is shown in the right panel, where we report the relative error between the 
code results and the simple analytical relation I I78II . In both panels, the horizontal axis displays the confined nucleons fraction Xc at 
r = Rcc- A horizontal dotted line marks the 2% error threshold in the right panel. 



between the free superfluid neutrons and the crust is restricted. Therefore, the shear mode frequencies tend to the values 
expected for models without a superfluid component. This effect is clear when we consider the plane-wave approximation (I70II 
in the limit ^ m. In this case, we have e* ~ Xc and equation (|70p becomes (jAndersson et al.ii2009i ) 

(77) 



which is the formula expected for pure elastic crust without superfluid neutrons (iMcDermott et al.lll988D . For the stellar 
models and effective masses considered in this work, there is roughly a 10% difference between the shear mode spectrum 
determined with and without a superfluid component in the crust. 

Finally, we study the effects of composition stratification on the spectrum. We consider a sequence of models D, where 
we set /t* = 10~* and take the entrainment to vanish. The original plane-wave relation (|7Up is (obviously) not adapted to 
describe the shear modes in stratified stars. However, we obtain a useful description of the spectrum from equation (|75p . 
according to which we have; 

a{xc) = ^f^, (78) 

\ \Xc)av 

where ai is the frequency for the non-stratified, {xc)av = 1, case. Figure |6] shows the frequencies of the first six shear 
modes for the set of D models. In the left panel, the numerical eigenfrequencies are compared to the values estimated from 
equation (|78p . The accuracy of the averaged plane- wave relation decreases when Xc becomes smaller at the crust/core interface. 
This behaviour is better illustrated in the right panel of Fig. [6l where we show the relative error of the approximation. In 
particular, we note that equation (|75p provides more accurate results for the realistic entrainment configuration effects than 
for the composition gradients. Therefore, when both the component fraction and the entrainment vary, the shear modes can 
be described by equation (|75p with accuracy limited by the composition gradient effects. 

Tables [2] and [3] provide the eigenfrequencies of the I = 2 interface mode and the first three I = 2 spheroidal shear modes 
for some selected models A and D (see the table captions for details). In Tabled the mode frequencies are given for different 
confined nucleon fractions Xc- In Table [S] the stellar models include the realistic entrainment profile. In order to test these 
results, we can convert the dimensionless eigenfrequencies of Table [2] into physical units and compare them to results from 
literature. As far as we know, there are no published results for shear and interface modes in single fluid stars with polytropic 
EoS and Newtonian gravity. However, there have been a number of analyses of the shear and interface mode spectrum for 
various other neutron star models. We can asses to what extent our results are consistent with these studies. To this end, we 
select from the stellar models A and D described in Sec. [3] the model with the results in the literature rely on the 

single fluid "approximation". From Tableland the central mass densities calculated in Sec. [3l the oscillation frequencies of 
the ^si and ^S2 modes are, respectively, u — 534.68 Hz and u — 986.49 Hz for model A, and u = 444.69 Hz and u = 812.43 Hz 
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Table 2. Frequencies of the I = 2 interface mode and the first three I = 2 spheroidal and toroidal shear modes for some of the stellar 
models studied in this work. Frequencies are given in units of a/^Gpo, where G is the gravitational constant and po represents the 
central mass density. The first column labels the stellar model. The non-stratified models A have constant proton fraction, Xp = 0.1 and 
the confined nucleon fraction Xc shown in the second column. The stratified models D have the same central proton fraction, Xp = 0.2 
at r = 0, while the values of Xc at r = Rcc are given in the second column. All models have zero entrainment and /t* = 10~^. 



Model 




2i 


2si 


2S2 




Ho 


^ti 


2t2 


A 


0.1 


0.1195 


1.2008 


2.2128 


3.2143 


0.0677 


1.1977 


2.2041 


A 


0.5 


0.0585 


0.4248 


0.7824 


1.4377 


0.0303 


0.5356 


0.9857 


A 


1.0 


0.0440 


0.3793 


0.6998 


1.0165 


0.0215 


0.3807 


0.7005 


D 


0.1 


0.0796 


0.5858 


1.0160 


1.4539 


0.0525 


0.5846 


1.0138 


D 


0.5 


0.0589 


0.3729 


0.6722 


0.9697 


0.0271 


0.3723 


0.6718 


D 


1.0 


0.0607 


0.3236 


0.5912 


0.8537 


0.0219 


0.3256 


0.5911 



Table 3. This table provides the frequencies of the I = 2 inte rface mode and the first three spheroidal and toroidal shear modes for 
models A and D with the "realistic" entrainment determined bv lCarter et al. I l l2005t ): [Ch"ame]| | |2005| . [2006h . See Sec. l3.2l for more details. 
Frequencies are given in dimensionless units, a/y/Gpo- The two stellar configurations selected in this table are: i) a model A with constant 
Xp = 0.1 and Xc = 0.3, and ii) a model D with Xp = 0.2 at r = and Xc = 0.3 at r = Rcc- All models have /t* = 10~*. 



Model 


2i 


2si 


2S2 




ho 


2ti 




A 


0.0431 


0.4265 


0.7765 


1.1227 


0.0391 


0.6915 


1.2726 


D 


0.0422 


0.3793 


0.6305 


0.8941 


0.0238 


0.3556 


0.6226 


for model D. For the I = 2 interface 


mode, 


we find 


= 62.03 


Hz and v 


= 83.41 


Hz for 


models A 



mode frequencies are consistent witli the results obtained by iMcDermott et al 



5.1.2 Torsional modes 

Motivated by ma^netar seismology, the spect rum of toroidal shear modes has already been studied in relativistic two-fluid 



neutron stars i Samuelsson fc Andersson 20091 ) . Our Newtonian analysis may therefore seem somewhat superfluous. However, 



if we want to improve our current stellar models by considering the effects of rotation and magnetic field, the perturbation 
problem becomes more complex (as the dynamics couples the polar and the axial perturbations). As the problem is more 
tractable in the context of Newtonian gravity, it makes sense to provide the corresponding non-rotating results here. 

Following the strategy adopted for the spheroidal shear modes, we study the dependence of the I = 2 toroidal modes on 
the shear modulus, composition fraction and entrainment. As previously, the plane-wave relation (|70p describes the spectrum 
of the non-stratified models A accurately, while the averaged version 1)75^ is more appropriate for both models D and models 
with realistic entrainment. It is not surprising that the main conclusions drawn for the s-modes remain valid also for the 
t-modes, as both sets of oscillations are restored by the crust elasticity. Given the similarities to the s-modes, we do not 
illustrate the torsional mode properties further, but the frequencies of the first three / = 2 t-modes are given in Tables [5] 

and [31 

The frequency of the fundamental torsional mo de, "^to, is typically between 10 and 50 Hz ( McDermott et al. 19881 : 



Samuelsson fc Andersson 20091 : Steiner fc WattsI 20091 ) . depending on the chosen model for the crust. Using the stellar pa- 



rameters in Table [2] we find that v = 30.31 Hz for model A and 30.09 Hz for model D. These values depend on the shear 
modulus as jxV^ , and hence the frequencies assume lower values for /i* < 10~^. The results are certainly consistent with the 
frequencies expected for the fundamental torsional modes. 



5.2 Acoustic modes 



Superfluid two-constituent stars can sustain two classes of fundamental and pressure modes, the ordinary and superfluid 
modes, respectively. The ordinary modes, for which the two fluid components tend to move in phase, have spectral properties 
similar to the f- and p-modes of single fluid stars. Given previous results fo r that problem, we expect their frequencies to be 
weakly affected by the presence of the elastic crust ( McDermott et al. 19881 ). The effects of the crust on the superfluid modes. 



for which the fluid elements oscillate in counter-phase, have not been studied previously. In this section, we provide the first 
analysis of this problem. We begin by considering the acoustic spectrum of non-stratified stellar models and then extend the 
analysis to stars with composition gradients. 
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Figure 7. Ordinary and superfluid 1 = 2 acoustic modes for models A with zero entrainment, Xp = 0.1 and /t* = lO"**. The left panel 
displays the dependence of the f- and p-mode frequencies on the crust composition Xc for a model A with Rcc = 0.9iJ. The right panel 
shows the superfluid I = 2 {- and p-modes for models A with different Xc and crust thickness, i.e. crust/core interface position Rcc (see 
legend). The mode frequencies are shown in dimensionless units. For non-stratified models, the ordinary modes depend weakly on the 
crust composition (left panel), while the superfluid mode frequencies change significantly with respect to Xc and Rcc- 



Table 4. Frequencies of the ordinary and superfluid 1 = 2 acoustic modes for a selection of stellar models A and D. Frequencies are 
given in units of rr/^Gpo, where G is the gravitational constant and po represents the central mass density. The flrst column labels the 
stellar model. The non-stratified models A have constant proton fraction, Xp = 0.1 and the confined nucleon fraction Xc is shown in the 
second column. The stratified models D have the same central proton fraction, Xp = 0.2 at r = 0, while the value of Zc at r = Rcc is 
given in the second column. All models have zero entrainment and /t* = 10~^. 



Model 


Xc 


2£0 


2pO 


2„o 
P2 


2_,o 
P3 


2£S 


2p| 


2_s 
P2 


2_s 

P3 


A 


0.1 


1.3848 


3.9069 


6.1122 


8.1952 


1.8687 


4.2831 


6.4227 


8.4708 


A 


0.5 


1.3848 


3.9071 


6.1123 


8.1949 


2.4209 


4.2727 


6.1482 


9.0163 


A 


1.0 


1.3848 


3.9073 


6.1125 


8.1948 


2.5398 


4.2689 


6.0909 


9.2041 


D 


0.1 


1.3085 


3.5649 


6.4524 


8.4873 


1.9319 


4.3871 


6.8852 


9.2062 


D 


0.5 


1.3462 


4.3439 


7.2824 


8.1749 


2.2839 


3.4936 


6.1328 


9.1061 


D 


1.0 


1.3392 


3.6177 


7.8464 


10.2406 


2.3946 


5.2819 


6.1476 


8.6563 



For a non-stratified model A with standard parameters, i.e. zero entrainment, Xp = 0.1, Rcc = 0.97? and p,* = 10~ , 
we determine tlie I = 2 i- and p-modes for different confined nucleon fractions in tlie crust, Xc. Results for the ordinary and 
superfluid modes are shown in the left panel of Fig. [7| and in Tabled] The dependence of the ■^f° and '^p" modes on Xc is very 
weak, as expected. For instance, the variation of the p-mode frequency for the two models, respectively, with Xc =0.1 and 
Xc = 1.0 is less than 0.01%. Focusing our attention on the ordinary modes, we also study the effects of the crust thickness on 
the eigenfrequencies. We consider a model A with Xp = 0.1 and Xc = 0.2, and change the interface location from Rcc = O.IR 
to Rcc ~ 0.99R. The deviation away from the ordinary mode frequency is very small even in the extreme cases, less than 
0.05%. 

The behaviour of the superfluid modes is different. Figure [7] and Table |3] show that the ^f and ^p^ modes depend 
signiflcantly on Xc, although in different ways. For instance, the frequencies of the ^f° and ^pj (^P2) modes are increasing 
(decreasing) functions of Xc, while the frequency of the ^p^ mode remains almost constant. The behaviour becomes even more 
complex when we vary the crust thickness. In the right panel of Fig. [71 we show the superfluid I = 2 f- and p-modes for three 
models A with different crust thickness, i.e. Rcc/R = 0.85,0.90 and 0.95, respectively. The crust region clearly affects the 
spectrum of superfluid modes and it is not (yet) clear to us whether a simple relation can be found to describe the observed 
behaviour. For tiie ^f' mode, we see that thicker crusts and larger fractions of confined nucei both tend to increase the mode 
frequency. 

In order to understand the ^f° mode better, we study its dependence on the change in composition (due to the chosen 
chemical gauge) between the core and crust, expressed in terms of A For a given Rcc, we construct two sets of 

models A. In the first set, we choose Xp = 0.1 and vary Xc, while in the second model we choose Xc = 0.9 and vary Xp. We 
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Figure 8. This figure displays tlie dependence of the superfluid I = 2 f-mode frequency on the crust/interface position Rcc and on the 
composition step Ax = Xc — Xp. The dimensionless eigenfrequencies are shown on the vertical axis. The background star is a model A 
with zero entrainment and /i* = 10~^. In the left panel, the -^f mode frequencies are determined for two different sequences of models 
A and for four different crust thickness, Rcc (see legend). The first sequence of stellar models has constant Xp = Xp = 0.1 and different 
Xc, while the second sequence has constant Xc = 1 — a;* = 0.9 and different Xp. The quantity x* is a parameter that identifies the two 
sequence of models A. The mode frequencies for the first sequence are denoted with dashed-lines, while those for the second sequence 
with symbols. The right panel shows the effect of Ax and Rcc on the mode frequencies for two sets of models A with x* = 0.5. This 
figure shows that for any two families of models A, denoted by x*, the eigenfrequencies are described by the same function of A, i.e. 
O- = cTaip (Ax). 



find that for both these classes of models, the mode frequency is described by the same function of Aa;. Results for this 
function, which we denote a = a (Ax), are shown in the left panel of Fig. |8] 

We find a similar behaviour if we extend this analysis to other stellar models, e.g., a sequence of models A with constant 
Xp = Xp and a second sequence with constant Xc — 1 ~ Xp. The parameter a;* identifies the two sets of models. For a given 
Xp, the mode frequency has the same functional dependence on Aa;. Note that the frequencies are different for distinct 
classes of stellar models, i.e. depend also on the parameter Xp, 

a = (Aa;) . (79) 

In the right panel of Fig.[8l we show the mode spectrum for two families of models A with x* = 0.5. In this case, we vary the 
composition fraction and the crust thickness extending the parameter space towards unrealistic configurations {Rcc ~ 0-lR) 
in order to highlight the qualitative behaviour. 

Let us now consider the acoustic mode spectrum of stratified stellar models. For models D with zero entrainment, we 
construct a sequence of stellar configurations with different composition Xc at the crust /core interface. The frequencies of the 
I = 2 f-mode and the first three p-modes for these models are shown in Fig. |9] and Table [l] The main effect of composition 
stratification is the coupling between the co- and counter-moving degrees of freedom of the oscillations. The interaction between 
ordinary and superfluid modes is evident in the spectrum (see Fig. [9]), where avoided crossings appear as the composition 
is varied. The chemical coupling also modifies the frequencies of the ordinary modes, as seen from the data in Table (4] For 
the ^f° mode, there is a frequency variation of 2.8% between models D with Xc — 0.1 and Xc — 1.0. However, this change 
may be due to the properties of the D models being somewhat different. In order to impose the continuity conditions at the 
crust/core boundary, the total mass density of the crust is a function of Xc (see Table [T|. This mass variation of the crust 
can modify the ^f° mode, as its frequency is proportional to the mean density of the stellar model, ~ M/R^ . We have not 
investigated this further, since the determined changes in the mode frequencies is not very significant. 

Now we turn to the effects that the large effective mass expected in the inner crust has on the acoustic spectrum. We 
consider a set of models D with different composition Xc at the crust /core interface and determine the frequencies of the 
ordinary and superfluid I = 2 fund amental mode s. In the right panel of Fig. [9] we compare the mode frequencies for the 
vanishing and realistic entrainment (|Chamelll2005l ) cases. The results show that the effective mass has a noticeable impact 
on the frequency of the fundamental mode. For the mode the frequency changes signiflcantly for the entire sequence of 
models D, while for the ^f ° mode we flnd a variation of 12.4% for a model D with Xc = 0.7, but only 0.01% for Xc — 0.3. This 
is an interesting result that deserves further investigation. However, in order to study more realistic neutron star models we 



need t o work in the framework of General Relativity. Developments in this direction are in progress, building on iLin et al. 



I 2008ll . 
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Figure 9. Frequencies of ttie ordinary and superfluid I = 2 acoustic modes for stellar models with composition stratification. The left 
panel shows the mode frequency dependence on the crust component Xc determined at the crust /core interface, r = Rcc- On the vertical 
axis the frequencies are in dimensionless units. The background model D has zero entrainment, /i* = 10"'* and Rcc = 0.88R. The p-mode 
spectrum exhibits avoided crossings between the superfluid and the ordinary modes. These arise from the interaction between the co- 
and counter-moving degrees of freedom that are coupled in stellar models by composition gradients. The right panel shows the effect 
of a "realistic" entrainment on the f-mode frequencies. The solid lines represent the results for zero entrainment models (the same as 
in the left panel), while the dashed lines denote the frequencies for models D with the Chamel entrainment configuration. The large 
entrainment of the deep regions of the crust affect both the superfluid and the ordinary 1 = 2 f-mode. We find that the variation of the 
■^f° mode frequency is about 12.4% for a model D with Xc = 0.7, but decreases to 0.01% for Xc = 0.3 



F inally, we test the result s against the acoustic spectrum o btained from time evolutions of superfluid neutron stars without 
crust /Pass amonti et al 20091 : Passamonti &: Andersson 2011 ). For the non-rotating model A, the frequencies of the ordinary 
and superfluid acoustic modes agree to better than 0.5%. This comparison is limited by the accuracy of the mode-extraction 
from the time-evolutions, so the result is acceptable. 



6 CONCLUSIONS 

In this paper we have studied the efi'ects of an elastic crust on the oscillation spectrum of superfluid neutron stars. Buildin g 



on recent results for the perturbations of two-fluid systems with both elasticity and superfluidity l| Andersson et al.l 120111 ), 
we considered Newtonian models of mature neutron stars. The stellar core was modelled as a mixture of superfluid neutrons 
and a conglomerate of charged particles, while the inner crust was described by a lattice of nuclei permeated by a gas of 
free superfluid neutrons. Focussing on nonrotating superfluid stars, we studied for the first time the effects of elasticity, 
entrainment and composition stratification on both shear modes and acoustic modes. Our results show that the superfluid 
neutrons in the crust may have considerable impact on the star's spectrum. The superfluid component clearly cannot be 
neglected if we want to determine an accurate oscillation spectrum of mature neutron stars. If we compare stellar models with 
and without superfluid neutrons in the crust, the difference between the shear mode frequencies is very large, but becomes 
smaller (about 10%) when we consider a "realistic" effective neutron mass. A superfluid component in the crust affects also the 
frequency of the superfluid acoustic modes, which are mainly dominated by a counter-phase motion between the two matter 
constituents. These superfluid modes are very sensitive to the composition difference between the core and the crust, and to 
the crust thickness. Ordinary acoustic modes, in which the stellar constituents oscillate mainly in phase, are less affected by the 
presence of the crust. However, the spectrum of ordinary modes becomes more complex when we consider stellar models with 
composition stratiflcation and more "realistic" entrainment. The spectrum exhibits avoided crossings between superfluid and 
ordinary pressure modes due to the chemical interaction between the co- and counter-moving degrees of freedom. Furthermore, 
for some specific stellar models we have found that a large effective mass may affect also the frequency of the ordinary f-mode 
at the 10% level. 

The results presented in this paper are important as they represent the first analysis of the polar modes of the coupled 
elastic superfluid system, complementing previous work on the corresponding axial mode problem. Our understanding of the 
latter problem remains more detailed, as it has also been considered in General Relativity. This is important since one must 
work in the context of relativity in order to be able to consider realistic neutron star equations of state in a meaningful way. The 
results we have presented were obtained for "relatively" simple, more phenomenological, stellar models. These models have the 
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advantage that one can explore the dependence on the key parameters, e.g., the composition, but they are obviously not "real- 
istic" . Nevertheless, the insights gained from the analysis provide useful guidelines for developing and testing fully relativistic 
models in the future. The state-of-play as far as such developments are concerned is tha t the required linear perturbation 
framework is complete and ready to be used in various applications of astrophysical interest ( Andersson. Samuelsson fc HaskeH 



2011f ). This step is, obviously, non-trivial but we expect to make progress in the near future. In order to make these models 



truly realistic, there also needs to be progress on various equation of state issues. We need more detailed models for the 
superfluid entrainment and the related effective masses, especially in the crust region. Further work in this direction should 
certainly be encouraged. After all, the entrainment is one of the key parameters that we may be able to infer from future 
neutron star seismology. 
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APPENDIX A: PERTURBATION EQUATION 

In this section, we provide the complete linearised equations for the polar sector. 



Al Core 

In the core, the eigenvalue problem can be studied by the set of perturbation variables defined in equation (|4ip. which we 
re-state here for completeness: 
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These dimensionless functions obey a system of six linearised equations: 
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In the crust, we use instead the dimensionless variables defined in equation (|43p . i.e., 
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where in this work qc — po- 
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The system of equations that describe the oscillation modes of a two fluid system in the crust is given by 
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where the have defined the following quantities: 
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In equation (|A14[) . the perturbation of the total mass density has the following form: 
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